/************************************************************************
Copyright (c) 2010, 2013, C.R. All rights reserved. 

E-mail: crtoyou@gmail.com; crtoyou@163.com

You can redistribute it and/or modify it under the terms of the
GNU Lesser General Public License as published by the
Free Software Foundation version 3 of the License.

This source code is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. if not, write to the Free Software Foundation, Inc., 
51 Franklin St, Fifth Floor, Boston, MA 02110-1301  USA
************************************************************************/

#ifndef __CR_MATRIX_WITH_DET_H__
#define __CR_MATRIX_WITH_DET_H__

#include <crbsplus.h>

namespace cr {

	template <class ele_type, short ROW, short COL = ROW>
	class matrix
	{
	public:
		const short cnt_row, cnt_col;

	public:
		matrix() : cnt_row(ROW), cnt_col(COL)
		{
			ASSERT(cnt_row>0 && cnt_col>0);
			for (int i=0; i<ROW; ++i)
				for (int j=0; j<COL; ++j)
					load[i][j] = 0;
		}

		matrix(const ele_type arr[ROW][COL])
			: cnt_row(ROW), cnt_col(COL)
		{
			ASSERT(cnt_row>0 && cnt_col>0);
			for (int i=0; i<ROW; ++i)
				for (int j=0; j<COL; ++j)
					load[i][j] = arr[i][j];
		}

		matrix(const matrix<ele_type, ROW, COL>& other)
			: cnt_row(ROW), cnt_col(COL)
		{
			ASSERT(cnt_row>0 && cnt_col>0);
			for (int i=0; i<ROW; ++i)
				for (int j=0; j<COL; ++j)
					load[i][j] = (ele_type)other[i][j];
		}

		matrix<ele_type,ROW,COL>& operator=(const matrix<ele_type,ROW,COL>& other)
		{
			ASSERT(cnt_row>0 && cnt_col>0);
			for (int i=0; i<ROW; ++i)
				for (int j=0; j<COL; ++j)
					load[i][j] = (ele_type)other[i][j];
			return *this;
		}

	protected:
		typedef ele_type row_type[COL];

		inline row_type& getrow(short r)
		{
			ASSERT(r>=0 && r<cnt_row);
			return load[r];
		}

		inline const row_type& getrow(short r) const
		{
			ASSERT(r>=0 && r<cnt_row);
			return load[r];
		}

		inline ele_type& getele(row_type& row, short c)
		{
			ASSERT(c>=0 && c<cnt_col);
			return row[c];
		}

		inline const ele_type& getele(const row_type& row, short c) const 
		{
			ASSERT(c>=0 && c<cnt_col);
			return row[c];
		}

		class row_vector 
		{
		public:
			row_vector(matrix<ele_type, ROW, COL>& mat, short r)
				: mat_ref(mat), r_ref(r)
			{}

			inline ele_type& operator[](short c)
			{
				return mat_ref.getele(mat_ref.getrow(r_ref), c);
			}

			inline const ele_type& operator[](short c) const
			{
				return mat_ref.getele(mat_ref.getrow(r_ref), c); 
			}

		private:
			matrix<ele_type, ROW, COL>& mat_ref;
			const short r_ref;
		};

	public:
		inline row_vector operator[](short r)
		{
			return row_vector(*this, r);
		}

		inline const row_vector operator[](short r) const
		{
			return row_vector(*const_cast<matrix<ele_type,ROW,COL>*>(this), r);
		}

	protected:
		ele_type load[ROW][COL];


		//////////////////////////////////////////////////////////////////////////
		//
		//nested class submatrix declaration:
		//
		//construct submatrix out of matrix is forbidden 
		//since the referred matrix 
		//may but forbid to dis-construct before its submatrix
		//////////////////////////////////////////////////////////////////////////

	protected:
		// !! the ROW and COL are for main matrix. 
		template <class ele_type, short ROW, short COL = ROW> 
		class submatrix 
		{
		public:
			const short cnt_row, cnt_col;

		public:
			submatrix(const matrix<ele_type,ROW,COL>& owner)
				: mat(owner)
				, mask_v(true), mask_h(true)
				, cnt_row(ROW), cnt_col(COL)
			{
				ASSERT(cnt_row>0 && cnt_col>0);
			}

			submatrix(const submatrix& parent, short i_row, short  i_col)
				: mat(parent.mat)
				, mask_v(parent.mask_v, i_row), mask_h(parent.mask_h, i_col)
				, cnt_row(parent.cnt_row-1), cnt_col(parent.cnt_col-1)
			{
				ASSERT(cnt_row>0 && cnt_col>0);
			}

			submatrix(const submatrix& other)
				: mat(other.mat)
				, mask_v(other.mask_v), mask_h(other.mask_h)
				, cnt_row(other.cnt_row), cnt_col(other.cnt_col)
			{
				ASSERT(cnt_row>0 && cnt_col>0);
			}

			//submatrix& operator=(const submatrix& other)
			//	: mat(other.mat)
			//	, mask_v(other.mask_v), mask_h(other.mask_h)
			//	, cnt_row(other.cnt_row), cnt_col(other.cnt_col)
			//{
			//	ASSERT(cnt_row>0 && cnt_col>0);
			//	return *this;
			//}

		protected:
			typedef ele_type row_type[COL];

			inline const row_type& getrow(short r) const
			{
				ASSERT(r>=0 && r<cnt_row);
				return mat.getrow(mask_v.seek(r));
			}

			inline const ele_type& getele(const row_type& row, short c) const
			{
				ASSERT(c>=0 && c<cnt_col);
				return mat.getele(row, mask_h.seek(c));
			}

			class subrow_vector
			{
			public:
				subrow_vector(const submatrix<ele_type, ROW, COL>& sub, short r)
					: sub_ref(sub), r_ref(r)
				{}

				inline const ele_type& operator[](short c) const
				{
					return sub_ref.getele(sub_ref.getrow(r_ref), c); 
				}

			private:
				const submatrix<ele_type, ROW, COL>& sub_ref;
				const short r_ref;
			};

		public:
			inline const subrow_vector operator[](short r) const 
			{
				return subrow_vector(*this, r);
			}

			template <short NEW_ROW, short NEW_COL>
			operator matrix<ele_type, NEW_ROW, NEW_COL>()
			{
				matrix<ele_type, NEW_ROW, NEW_COL> mat;
				ASSERT(cnt_row == mat.cnt_row);
				ASSERT(cnt_col == mat.cnt_col);
				for (short r=0; r<cnt_row; ++r)
					for (short c=0; c<cnt_col; ++c)
						mat[r][c] = (*this)[r][c];
				return mat;
			}

		private:
			friend class matrix<ele_type, ROW, COL>;

			//////////////////////////////////////////////////////////////////////////
			// all the elements are reference to its owner. 
			// co gradient with its owner. 
			//////////////////////////////////////////////////////////////////////////

			const matrix<ele_type,ROW,COL>& mat;
			const bitmask<ROW> mask_v;
			const bitmask<COL> mask_h;
		};


	public:
		template <short OLD_ROW, short OLD_COL>
		matrix(const submatrix<ele_type, OLD_ROW, OLD_COL>& sub)
			: cnt_row(ROW), cnt_col(COL)
		{
			ASSERT(cnt_row == sub.cnt_row);
			ASSERT(cnt_col == sub.cnt_col);
			for (short r=0; r<cnt_row; ++r)
				for (short c=0; c<cnt_col; ++c)
					load[r][c] = sub[r][c];
		}

	};


	template <class ele_type, short ROW, short COL = ROW> 
	class matrixmaths : public matrix<ele_type, ROW, COL>
	{
	public:
		matrixmaths() : matrix()
		{}

		matrixmaths(const ele_type arr[ROW][COL])
			: matrix(arr)
		{}

	public:
		bool det(ele_type* val) const 
		{
			ASSERT(val);
			submatrixmaths<ele_type,ROW,COL> sub(*this);
			return sub.det(val);
		}

	protected:
		template <class ele_type, short ROW, short COL = ROW> 
		class submatrixmaths : public submatrix<ele_type, ROW, COL>
		{
		public:
			submatrixmaths(const matrixmaths<ele_type,ROW,COL>& owner)
				: submatrix(owner)
			{}

			submatrixmaths(const submatrixmaths& parent, short i_row, short  i_col)
				: submatrix(parent, i_row, i_col)
			{}

		public:

			//////////////////////////////////////////////////////////////////////////
			//!scan only and any row or column is ok. 
			//it is preferred to choose the last row to scan. 
			//scan last row makes the masks are begin from first bit, 
			//so that promote the performance of masks. 
			//////////////////////////////////////////////////////////////////////////

			bool det(ele_type* val) const 
			{
				const short& cnt = cnt_row;
				ASSERT(val);
				ASSERT(cnt_col == cnt_row); //only spare matrix has det. 
				ASSERT(cnt > 0);
				ASSERT(cnt < 20);//unexpected 
				*val = 0;
				if (1 == cnt) {
					*val = (*this)[0][0];
				}
				else if (2 == cnt) {
					*val = (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0];
				}
				else if (3 == cnt) {
					*val += (*this)[0][0] * (*this)[1][1] * (*this)[2][2]; //++plus++
					*val += (*this)[0][1] * (*this)[1][2] * (*this)[2][0]; //++plus++
					*val += (*this)[0][2] * (*this)[1][0] * (*this)[2][1]; //++plus++
					*val -= (*this)[0][2] * (*this)[1][1] * (*this)[2][0]; //--minus--
					*val -= (*this)[0][1] * (*this)[1][0] * (*this)[2][2]; //--minus--
					*val -= (*this)[0][0] * (*this)[1][2] * (*this)[2][1]; //--minus--
				}
				else {
					//scan last row, walk column from 0 to cnt_col. 
					short r = cnt - 1;
					for (int c=0; c<cnt; ++c) {
						submatrixmaths<ele_type, ROW, COL> sub(*this, r, c);
						ele_type subval;
						if (!sub.det(&subval))
							return false;
						if ((r+c)%2) {
							*val -= (*this)[r][c] * subval;
						}
						else {
							*val += (*this)[r][c] * subval;
						}
					}
				}
				return true;
			}
		};

	};

	typedef matrix<int, 4, 4> matrix4_d;
	typedef matrix<int, 3, 3> matrix3_d;
	typedef matrix<int, 2, 2> matrix2_d;

	typedef matrixmaths<int, 4, 4> matrixmaths4_d;
	typedef matrixmaths<int, 3, 3> matrixmaths3_d;
	typedef matrixmaths<int, 2, 2> matrixmaths2_d;

	typedef matrix<double, 4, 4> matrix4_f;
	typedef matrix<double, 3, 3> matrix3_f;
	typedef matrix<double, 2, 2> matrix2_f;

	typedef matrixmaths<double, 4, 4> matrixmaths4_f;
	typedef matrixmaths<double, 3, 3> matrixmaths3_f;
	typedef matrixmaths<double, 2, 2> matrixmaths2_f;

}//namespace

#endif //__CR_MATRIX_WITH_DET_H__